Multilevel model
Jiayue Ao
2023-07-10
1. Packages
# Packages
library(haven)
library(tidyverse)
library(naniar)
library(scales)
library(lme4)
library(jtools)
library(stargazer)
library(sjPlot)
library(performance)
library(nlme)
library(lattice)
library(sjmisc)
library(sjPlot)
library(glmmTMB)
library(psych)2. Dataset
2.1 Data Source
Data Source: Understanding Society – The UK Household Longitudinal Study, “w_indresp.dta”, w means wave a to wave l (expect wave h) Download from: UK Data Service-SN6614-Understanding Society: Waves 1-12, 2009-2021 and Harmonised BHPS: Waves 1-18, 1991-2009
2.2 Longitudinal dataset
The data are from the Understanding Society The UK Longitudinal study, where this dissertation selected the survey data from wave1 in 2009 to wave12 in 2021 indresp.dta, except for wave 8 survey. Because the main data of political Interest were omitted in wave 8 survey data. Therefore, there are total 11 waves data needed to study. The following most important thing to do is to combine the datasets from the 11 independent files into one longitudinal dataset.
Based on the definition of older people, the project will select the observations aged between 60 and 69 in the first wave and attend all 11 waves(Fletcher, 2021).
2.2.1 Selecting variables
The project will collect the datat from the Understanding Society “w_indresp.dta”. Because the data is scattered in 11 separate dta files, the project needs to be merged into one longitudinal dataset.
Based on the research question, the project will select 1 dependent variable and 5 independent variables. The dependent variable named political interest. In the survey, political interest is expressed as “_vote6”. The value 1 means the respondent does not interested in politics. The value 2 means the respondent not very interested in political interest. The value 3 means the respondent fairly interested in politics. The value 4 means the respondent very interested in politics. The first independent variable is age. In the survey, age is expressed as “_dvage”. The second independent variable is retirement. In the survey, retirement is expressed as “a_retdatey”(wave 1), “_jbendreas”(wave2-6), “_jbendreas6”(wave7-11). The value 1 of retirement means the older people is retired. The value 0 means the older people does not retire. The third independent variable is health. In the survey, health status is expressed as “_health”. Health is a binary variable, value 1 means the observation is long-standing illness or disability, and value 0 means the observations does not long-standing illness or disability. For the variavle for the married status in the survey is “_marstat_dv”. The fourth independent variable widow and divorce could be found in “_marstat_dv”. The fifthe variable income is continuous varibale.
# Empty tibble
temp_data <- tibble()
# 11 years(waves)
year_letters <- c("a", "b", "c", "d","e","f","g","i","j","k","l")
year_nums <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
# Finding the variables from 11 data files
for (i in 1:length(year_letters)) {
year_letter <- year_letters[i]
year_num <- year_nums[i]
if (file.exists(paste0(year_letter, "_indresp.dta"))) {
year_data <- read_dta(paste0(year_letter, "_indresp.dta"))
if(year_num==1){
# Finding the variables
selected_variables <- year_data %>%
rename(pidp = pidp,
PI = paste0(year_letter, "_vote6"),
age = paste0(year_letter, "_dvage"),
employment = paste0(year_letter, "_employ"),
sex = paste0(year_letter, "_sex_dv"),
health = paste0(year_letter, "_health"),
edu = paste0(year_letter, "_qfhigh_dv"),
marital_status = paste0(year_letter, "_marstat_dv"),
income = paste0(year_letter, "_fimnnet_dv"),
retirement = a_retdatey) %>%
select(pidp, PI, age, employment, sex, health, edu,marital_status, income, retirement) %>%
mutate(wave = year_num,
retirement = ifelse(retirement>=0,1,NA))
}
if(year_num %in% 2:6){
# Finding the variables
selected_variables <- year_data %>%
rename(pidp = pidp,
PI = paste0(year_letter, "_vote6"),
age = paste0(year_letter, "_dvage"),
employment = paste0(year_letter, "_employ"),
sex = paste0(year_letter, "_sex_dv"),
health = paste0(year_letter, "_health"),
edu = paste0(year_letter, "_qfhigh_dv"),
marital_status = paste0(year_letter, "_marstat_dv"),
income = paste0(year_letter, "_fimnnet_dv"),
retirement = paste0(year_letter, "_jbendreas")) %>%
select(pidp, PI, age, employment, sex, health, edu,marital_status, income, retirement) %>%
mutate(wave = year_num,
retirement = ifelse(retirement==6,1,NA))
}
if(year_num %in% 7:11){
# Finding the variables
selected_variables <- year_data %>%
rename(pidp = pidp,
PI = paste0(year_letter, "_vote6"),
age = paste0(year_letter, "_dvage"),
employment = paste0(year_letter, "_employ"),
sex = paste0(year_letter, "_sex_dv"),
health = paste0(year_letter, "_health"),
edu = paste0(year_letter, "_qfhigh_dv"),
marital_status = paste0(year_letter, "_marstat_dv"),
income = paste0(year_letter, "_fimnnet_dv"),
retirement = paste0(year_letter, "_jbendreas6")) %>%
select(pidp, PI, age, employment, sex, health, edu,marital_status, income, retirement) %>%
mutate(wave = year_num,
retirement = ifelse(retirement==1,1,NA))
}
# Combine data together
temp_data <- bind_rows(temp_data, selected_variables)
}
}# Deal with retirement
temp_data = temp_data %>%
arrange(pidp,wave) %>%
group_by(pidp) %>%
fill(retirement,.direction = 'down') %>%
ungroup() %>%
mutate(retirement = ifelse(is.na(retirement),0,retirement))# Define NA missing value
temp_data$PI[temp_data$PI %in% c(-10, -9, -8, -7, -2, -1)] <- NA
temp_data$edu[temp_data$edu %in% c(-9, -8)] <- NA
temp_data$sex[temp_data$sex == -9] <- NA
temp_data$employment[temp_data$employment %in% c(-9, -8, -2, -1)] <- NA
temp_data$health[temp_data$health %in% c(-9, -8, -2, -1)] <- NA
temp_data$marital_status[temp_data$marital_status%in% c(-9, 0)] <- NA
temp_data$income[temp_data$income %in% c(-9, -4)] <- NA
# Delete NA missing value
temp_data <- na.omit(temp_data)# PI order change
# According to the survey design(Understanding Society, 2023), the lower numeric result shows that people more interested in PI. However, this is a departure from normal logic. The person might think that a higher number represents a higher PI.
# 1 is or not at all interested?
# 2 is not very
# 3 is Fairly
# 4 is Very
temp_data <- temp_data %>%
mutate(PI = case_when(
PI == 1 ~ 4,
PI == 2 ~ 3,
PI == 3 ~ 2,
PI == 4 ~ 1
))
# Education
# The line of demarcation is EDUCATION, so 0 is not graduate, 1 is graduate.
temp_data <- temp_data %>%
mutate(edu = case_when(
edu == 1 ~ 1,
edu == 2 ~ 1,
edu == 3 ~ 1,
edu == 4 ~ 1,
edu == 5 ~ 1,
edu == 6 ~ 1,
edu == 7 ~ 1,
edu == 8 ~ 1,
edu == 9 ~ 1,
edu == 10 ~ 1,
edu == 11 ~ 1,
edu == 12 ~ 1,
edu == 13 ~ 1,
edu == 14 ~ 1,
edu == 15 ~ 1,
edu == 16 ~ 1,
edu == 96 ~ 0
))
# Marital_status- 3 is widow, 4 is divorce
temp_data <- temp_data %>%
mutate(marital_status = case_when(
marital_status == 1 ~ 0,
marital_status == 2 ~ 0,
marital_status == 3 ~ 1,
marital_status == 4 ~ 2,
marital_status == 5 ~ 0,
marital_status == 6 ~ 0
))
# Setting Log_income
threshold <- 1e-10
# There existing the negative value in the income
temp_data <- temp_data %>%
mutate(income = ifelse(income <= 0, threshold, income))
# Log
temp_data <- temp_data %>%
group_by(pidp) %>%
arrange(wave) %>%
mutate(log_income = log(income)) %>%
ungroup()
# Setting time Lag t-1
temp_data <- temp_data %>%
group_by(pidp) %>%
arrange(wave) %>%
mutate(lag_retirement_1 = lag(retirement)) %>%
ungroup()
temp_data <- temp_data %>%
group_by(pidp) %>%
arrange(wave) %>%
mutate(lag_health_1 = lag(health)) %>%
ungroup()
temp_data <- temp_data %>%
group_by(pidp) %>%
arrange(wave) %>%
mutate(lag_marital_status_1 = lag(marital_status)) %>%
ungroup()
# Setting time Lag t-2
temp_data <- temp_data %>%
group_by(pidp) %>%
arrange(wave) %>%
mutate(lag_retirement_2 = lag(retirement,2)) %>%
ungroup()
temp_data <- temp_data %>%
group_by(pidp) %>%
arrange(wave) %>%
mutate(lag_health_2 = lag(health,2)) %>%
ungroup()
temp_data <- temp_data %>%
group_by(pidp) %>%
arrange(wave) %>%
mutate(lag_marital_status_2 = lag(marital_status,2)) %>%
ungroup()
# Categorical variables
temp_data$sex <- as.factor(temp_data$sex)
temp_data$edu <- as.factor(temp_data$edu)
temp_data$employment <- as.factor(temp_data$employment)
temp_data$marital_status <- as.factor(temp_data$marital_status)
temp_data$retirement <- as.factor(temp_data$retirement)
temp_data$health <- as.factor(temp_data$health)
temp_data$lag_marital_status_1 <- as.factor(temp_data$lag_marital_status_1)
temp_data$lag_marital_status_2 <- as.factor(temp_data$lag_marital_status_2)
temp_data$lag_retirement_1 <- as.factor(temp_data$lag_retirement_1)
temp_data$lag_retirement_2 <- as.factor(temp_data$lag_retirement_2)
temp_data$lag_health_1 <- as.factor(temp_data$lag_health_1)
temp_data$lag_health_2 <- as.factor(temp_data$lag_health_2)# Select aging population (first wave aged 60-69)
Longitudinal_data <- temp_data %>%
arrange(pidp, wave) %>%
group_by(pidp) %>%
filter(n() == 11,
first(age) >= 60,
first(age) <= 69) %>%
ungroup()# Recode the value variable
data <- Longitudinal_data %>%
ungroup() %>%
select(pidp:lag_marital_status_2) %>%
mutate(health = zap_label(health),
health = as.factor(health),
age = zap_label(age),
age = as.numeric(age),
income = zap_label(income),
income = as.numeric(income),
employment = fct_recode(employment,
'yes' = '1',
'no' = '2'),
employment = fct_relevel(employment,'no','yes'),
sex = fct_recode(sex,
'male' = '1',
'female' = '2'),
sex = fct_drop(sex),
health = fct_recode(health,
'yes' = '1',
'no' = '2'),
health = fct_relevel(health,'no','yes'),
edu = fct_recode(edu,
'not graduate' = '0',
'graduate' = '1'),
marital_status = fct_recode(marital_status,
'none' = '0',
'widow' = '1',
'divorce' = '2'),
retired = fct_recode(retirement,
'no' = '0',
'yes' = '1'),
lag_retired_1 = fct_recode(lag_retirement_1,
'no' = '0',
'yes' = '1'),
lag_retired_2 = fct_recode(lag_retirement_2,
'no' = '0',
'yes' = '1'),
lag_health_1 = fct_recode(lag_health_1,
'yes' = '1',
'no' = '2'),
lag_health_2 = fct_recode(lag_health_2,
'yes' = '1',
'no' = '2'),
lag_marital_status_1 = fct_recode(lag_marital_status_1,
'none' = '0',
'widow' = '1',
'divorce' = '2'),
lag_marital_status_2 = fct_recode(lag_marital_status_2,
'none' = '0',
'widow' = '1',
'divorce' = '2')) %>%
select(-retirement, -lag_retirement_1, -lag_retirement_2)3. Visualization of variables
vis_miss(data)
### 3.1 Politting age and politcial interest
data %>%
mutate(PI = as.factor(PI)) %>%
count(PI) %>%
ggplot(aes(x = PI, y = n, fill = PI)) +
geom_col() +
scale_fill_brewer(palette = 'RdBu',
labels = c("Or not at all interested", "Not very interested", "Fairly interested", "Very interested")) +
labs(x = 'Political Interest',
y = 'Frequency',
title = 'The Frequency of Political Interest') +
theme_bw() 
data %>%
mutate(PI = as.factor(PI)) %>%
count(PI, wave) %>%
ggplot(aes(x = wave, y = n, fill = PI)) +
geom_bar(stat = 'identity', position = 'fill') +
scale_x_continuous(breaks = 1:11) +
scale_y_continuous(labels = percent) +
scale_fill_brewer(palette = 'RdBu',
labels = c("Or not at all interested", "Not very interested", "Fairly interested", "Very interested")) +
labs(x = 'Waves',
y = 'Pencent (%)',
title = 'The Percentage of Each Degree of Political Interest in 11 Waves') +
theme_bw()+
theme_classic() +
theme(plot.title = element_text(size = 12))
data %>%
mutate(PI = as.factor(PI)) %>%
ggplot(aes(x = PI,y = age)) +
geom_violin(mapping = aes(fill = PI),alpha = 0.5) +
geom_jitter(mapping = aes(color = PI),alpha = 0.3) +
scale_fill_brewer(palette = 'RdBu',
labels = c("Or not at all interested", "Not very interested", "Fairly interested", "Very interested")) +
scale_colour_brewer(palette = 'RdBu',
labels = c("Or not at all interested", "Not very interested", "Fairly interested", "Very interested")) +
labs(x = 'Political Interest',
y = 'Ages',
title = 'The Degree of Political Interest at Different Ages') +
theme_bw()
plot_grpfrq(data$PI, data$age, type = "violin")
The figures shows that older people of different ages are concentrated
in different levels of political interest.
3.2 Plotting time-varing variable retirement and political interest
data %>%
count(PI,retired) %>%
ggplot(aes(x = PI,y = n,fill = retired)) +
geom_bar(stat = 'identity',position = 'fill') +
scale_y_continuous(labels = percent) +
labs(x = 'Political Interest',
y = 'Pencent (%)',
fill = 'Retirement Status',
title = 'The Percent of Political Interest at Different Retirement Status') +
theme_bw() +
theme(plot.title = element_text(size = 12))
3.3 Plotting income and political interest
data %>%
mutate(PI = as.factor(PI)) %>%
ggplot(aes(x = PI,y = income)) +
geom_violin(mapping = aes(fill = PI),alpha = 0.5) +
geom_jitter(mapping = aes(color = PI),alpha = 0.1) +
scale_fill_brewer(palette = 'RdBu',
labels = c("Or not at all interested", "Not very interested", "Fairly interested", "Very interested")) +
scale_colour_brewer(palette = 'RdBu',
labels = c("Or not at all interested", "Not very interested", "Fairly interested", "Very interested")) +
scale_y_sqrt() +
labs(x = 'Political Interest',
y = 'Income',
title = 'The Degree of Political Interest at Different Incomes') +
theme_bw() 
3.4 Plotting health status and political interest
data %>%
count(PI,health) %>%
ggplot(aes(x = PI,y = n,fill = health)) +
geom_bar(stat = 'identity',position = 'fill') +
scale_y_continuous(labels = percent) +
labs(x = 'Political Interest',
y = 'Pencent (%)',
fill = 'Health',
title = 'The Percent of Political Interest at Different Healths') +
theme_bw() +
theme(plot.title = element_text(size = 12))
3.5 Plotting widow and political interest
data %>%
count(PI, marital_status) %>%
ggplot(aes(x = PI, y = n, fill = marital_status)) +
geom_bar(stat = 'identity', position = 'fill') +
scale_y_continuous(labels = percent) +
labs(x = 'Political Interest',
y = 'Percent (%)',
fill = 'marital_status',
title = 'The Percent of Political Interest at Different marital_status') +
theme_bw() +
theme(plot.title = element_text(size = 12))
### 3.6 Descriptive Statistics
Descriptive <- summary(data)print(Descriptive) pidp PI age employment sex
Min. :1.021e+09 Min. :1.000 Min. :60.0 no :5184 male :2915
1st Qu.:1.360e+09 1st Qu.:2.000 1st Qu.:66.0 yes:1031 female:3300
Median :1.429e+09 Median :3.000 Median :69.0
Mean :1.436e+09 Mean :2.682 Mean :69.1
3rd Qu.:1.564e+09 3rd Qu.:3.000 3rd Qu.:72.0
Max. :1.634e+09 Max. :4.000 Max. :81.0
health edu marital_status income wave
no :3169 not graduate:1952 none :4703 Min. : 0 Min. : 1
yes:3046 graduate :4263 widow : 779 1st Qu.: 866 1st Qu.: 3
divorce: 733 Median : 1315 Median : 6
Mean : 1640 Mean : 6
3rd Qu.: 1996 3rd Qu.: 9
Max. :33184 Max. :11
log_income lag_health_1 lag_marital_status_1 lag_health_2
Min. :-23.026 yes :2763 none :4295 yes :2508
1st Qu.: 6.764 no :2887 widow : 686 no :2577
Median : 7.182 NA's: 565 divorce: 669 NA's:1130
Mean : 6.966 NA's : 565
3rd Qu.: 7.599
Max. : 10.410
lag_marital_status_2 retired lag_retired_1 lag_retired_2
none :3882 no :2598 no :2362 no :2126
widow : 600 yes:3617 yes :3288 yes :2959
divorce: 603 NA's: 565 NA's:1130
NA's :1130
4. A baseline variance components model
m0 <- lmer(PI ~ (1 | pidp), data = data)
summary(m0, correlation = FALSE)Linear mixed model fit by REML ['lmerMod']
Formula: PI ~ (1 | pidp)
Data: data
REML criterion at convergence: 11101.8
Scaled residuals:
Min 1Q Median 3Q Max
-4.9091 -0.4770 0.0257 0.5371 4.8796
Random effects:
Groups Name Variance Std.Dev.
pidp (Intercept) 0.5508 0.7421
Residual 0.2614 0.5112
Number of obs: 6215, groups: pidp, 565
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.68238 0.03189 84.12The aim of the study is to examine how the dependent variable PI depends on different pidp (that is, individual ids). The random intercept of the model only considers different pidp (individual ids) without any fixed effects. The model means that you are estimating the variance of the average PI value for each pidp, as well as the overall average of the PI values. The hierarchy of longitudinal data is panel data, Level 1: occasions, Level 2: pidp (Bell & Fairvrother & Jones, 2018).
REML criterion at convergence：11101.8 Random effects： -pidp (Intercept): the Variance of the random intercept for pidp is 0.5508, Std. Dev. is 0.7421. -Residual：the Variance of the model residual is 0.2614, Std. Dev. is 0.5112. The variance of the residuals (0.2614) is smaller than the variance of the pidp (0.5508). Therefore, The variation between older people is larger than the variation between occasions within older people.
Fixed effects: For constant term, estimate value 2.68238, Std. Error 0.03189, t value 84.12. The result of fixed effect means that the overall mean of the PI is 2.68238.
5. A linear growth curve model
5.1 Age, fixed explanatory variables
A first step in modelling the between occasion within older people, or level 1, variation would be to fit a fixed linear trend. Age as fixed explanatory variables.
m1 <- lmer(PI ~ age + (1 | pidp), data = data)
summary(m1, correlation = FALSE)Linear mixed model fit by REML ['lmerMod']
Formula: PI ~ age + (1 | pidp)
Data: data
REML criterion at convergence: 11067.3
Scaled residuals:
Min 1Q Median 3Q Max
-4.8734 -0.4631 0.0242 0.5295 4.8812
Random effects:
Groups Name Variance Std.Dev.
pidp (Intercept) 0.5476 0.7400
Residual 0.2595 0.5094
Number of obs: 6215, groups: pidp, 565
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.847213 0.127852 14.448
age 0.012086 0.001792 6.744REML criterion at convergence: 11067.3. The REML value is smaller compared to the previous model m0 without fixed explanatory variables. The resulte means that the model with fixed explanatory variables fits better.
Random effects: -pidp (Intercept): The variance of the pidp (Intercept) decreased slightly, from 0.5508 to 0.5476. The decrease is slight. -Resdisual: The variance of the residuals also decreased slightly, from 0.2614 to 0.2595. The variance of the random effects decreased slightly. Because the decrease was samll, suggesting that age did not explain a great deal of the variation within pidp.
Fixed effects: The intercept decreases from 2.69534 to 1.847213. The result means that the predicted value of the PI is 1.847213 when age is 0. The estimate of age is 0.012086. The result means that when the unit of age increase, the predicted value of the PI increases by 0.012086 units on average. Moreover, the t-value is 6.744, which is usually considered statistically significant, meaning that age is a significant predictor variable.
Likelihood Ratio Test (LRT) is equal to 34.51423, means the m1 fit is better
-2*(logLik(m0)-logLik(m1))'log Lik.' 34.51423 (df=3)5.2 Random intercept and random slope model
The section 5.2 will make the coefficient of age random at level 2. The model includes random intercept and random slope.
m2 <- lmer(PI ~ age + (1 + age | pidp), data = data)
summary(m2, correlation = FALSE)Linear mixed model fit by REML ['lmerMod']
Formula: PI ~ age + (1 + age | pidp)
Data: data
REML criterion at convergence: 10919.4
Scaled residuals:
Min 1Q Median 3Q Max
-4.9274 -0.4545 0.0394 0.4875 5.0774
Random effects:
Groups Name Variance Std.Dev. Corr
pidp (Intercept) 7.05370 2.65588
age 0.00128 0.03578 -0.96
Residual 0.23926 0.48914
Number of obs: 6215, groups: pidp, 565
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.843579 0.163468 11.28
age 0.012221 0.002288 5.34
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 1.31049 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?In the m2, the random slope is age. The definition of random slope is allowing the associations between variables to vary across higher-level entities (Bell & Fairbrothers &Jones, 2018).
Random Effect: - The variance of pidp (Intercept) is 7.05370 with a standard deviation of 2.65588, which means that there is significant variation in the baseline PI value (the value when age is 0) between different pidp groups. - The variance of age is 0.00128 with a standard deviation of 0.03578. this means that there is also significant variation in the slope of the effect of age on PI between different pidp groups. - The value of Corr is -0.96，and the value close to -1. Thus, the age as random slope is insignificant. Because everyone gets same coefficient for age and gets same random effect. Moreover, model failed to converge, and the model needs to be simplified. The study of the project focus on the group of older people. As the samples all about older people, thus the coefficient may not much different. Therefore, the following model will drop age as random slope and only include pidp as random intercept. The conclusion is different from the Bell & Fairbrothers &Jones (2018) that within-between RE model needs to include random slopes as much as possible.
Fixed Effect: - Intercept 1.843579 is the estimated average value of the PI when age is zero. Std. Error: The value of Std. Error is 0.163468. T-value is 11.28 which is statistically significant. -age The variable age is a significant predictor of the dependent variable. The coefficient for age is 0.012221. For every one-unit increase in age, the predicted value of PI increases by 0.012221 units, holding other factors constant. The Std. Error is 0.002288. T-value is 5.34. The result shows that the relationship between age and DV is statistically significant.
Likelihood Ratio Test (LRT) is equal to 147.8674.
-2*(logLik(m1)-logLik(m2))'log Lik.' 147.8674 (df=4)5.3 Predicted slopes
To generate the predicted slopes for each older people type
# Select the first 50 unique pidp
selected_pidp = unique(data$pidp)[30:50]
# Filter the data
selected_data = data[data$pidp %in% selected_pidp,]
# Generate fitted scores for the selected samples
selected_predscore = fitted(m2)[data$pidp %in% selected_pidp]
# Make the plot
xyplot(selected_predscore ~ age, data = selected_data, groups = pidp, type = c("p", "l"), col = "blue")
Since there are too many pidp (individuals), some of the pidp are
selected and the predictor scores of the predicted political interest
score are plotted. According to the picture, most of the pidp will have
a predictor value that rises with age.
5.4 Repeated measures modelling of non-linear polynomial growth
Check for non-linear relationships
m3 <- lmer(PI ~ age + I(age^2) + (1 | pidp), data = data)
summary(m3)Linear mixed model fit by REML ['lmerMod']
Formula: PI ~ age + I(age^2) + (1 | pidp)
Data: data
REML criterion at convergence: 11066.1
Scaled residuals:
Min 1Q Median 3Q Max
-4.8999 -0.4739 0.0197 0.5239 4.8617
Random effects:
Groups Name Variance Std.Dev.
pidp (Intercept) 0.5475 0.7399
Residual 0.2588 0.5087
Number of obs: 6215, groups: pidp, 565
Fixed effects:
Estimate Std. Error t value
(Intercept) -3.7566362 1.4259996 -2.634
age 0.1744738 0.0411951 4.235
I(age^2) -0.0011716 0.0002969 -3.946
Correlation of Fixed Effects:
(Intr) age
age -0.999
I(age^2) 0.996 -0.999# Select the first 50 unique pidp
selected_pidp = unique(data$pidp)[30:50]
# Filter the data
selected_data = data[data$pidp %in% selected_pidp,]
# Generate fitted scores for the selected samples
selected_predscore_2 = fitted(m3)[data$pidp %in% selected_pidp]
# Make the plot
xyplot(selected_predscore_2 ~ age, data = selected_data, groups = pidp, type = c("p", "l"), col = "blue")
The pciture shows that with the age increase, the political interest
predicted score will increase then decrase. The line of predicted score
is nolinear.
The likelihood statistic value is negative -146.7237. The result shows that the quadratic term would not improve the model.
-2*(logLik(m2)-logLik(m3))'log Lik.' -146.7237 (df=6)6. Full Multilvel Model
# Time-lag for income
data <- data %>%
group_by(pidp) %>%
arrange(wave) %>%
mutate(lag_income_1 = lag(log_income)) %>%
ungroup()
data <- data %>%
group_by(pidp) %>%
arrange(wave) %>%
mutate(lag_income_2 = lag(log_income,2)) %>%
ungroup()6.1 Full variables without time-lag
gm1 <- lmer(PI ~ age + retired + health + marital_status + health + log_income + (1 |pidp), data = data)
summary(gm1, correlation = FALSE)Linear mixed model fit by REML ['lmerMod']
Formula: PI ~ age + retired + health + marital_status + health + log_income +
(1 | pidp)
Data: data
REML criterion at convergence: 11086.4
Scaled residuals:
Min 1Q Median 3Q Max
-4.8791 -0.4674 0.0225 0.5293 4.8221
Random effects:
Groups Name Variance Std.Dev.
pidp (Intercept) 0.5401 0.7349
Residual 0.2596 0.5095
Number of obs: 6215, groups: pidp, 565
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.811367 0.133620 13.556
age 0.012884 0.001821 7.075
retiredyes -0.010130 0.063366 -0.160
healthyes -0.004139 0.017491 -0.237
marital_statuswidow -0.127886 0.043191 -2.961
marital_statusdivorce -0.050341 0.053953 -0.933
log_income 0.001519 0.002847 0.534summ(gm1,confint = T)MODEL INFO:
Observations: 6215
Dependent Variable: PI
Type: Mixed effects linear regression
MODEL FIT:
AIC = 11104.36, BIC = 11164.97
Pseudo-R² (fixed effects) = 0.01
Pseudo-R² (total) = 0.68
FIXED EFFECTS:
------------------------------------------------------------
Est. 2.5% 97.5% t val.
--------------------------- ------- ------- ------- --------
(Intercept) 1.81 1.55 2.07 13.56
age 0.01 0.01 0.02 7.07
retiredyes -0.01 -0.13 0.11 -0.16
healthyes -0.00 -0.04 0.03 -0.24
marital_statuswidow -0.13 -0.21 -0.04 -2.96
marital_statusdivorce -0.05 -0.16 0.06 -0.93
log_income 0.00 -0.00 0.01 0.53
------------------------------------------------------------
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
pidp (Intercept) 0.73
Residual 0.51
------------------------------------
Grouping variables:
-------------------------
Group # groups ICC
------- ---------- ------
pidp 565 0.68
-------------------------6.2 Time lag t-1
gm2 <- lmer(PI ~ age + lag_retired_1 + lag_health_1 + lag_marital_status_1 + lag_income_1 + (1 |pidp), data = data, na.action=na.exclude)
summary(gm2, correlation = FALSE)Linear mixed model fit by REML ['lmerMod']
Formula: PI ~ age + lag_retired_1 + lag_health_1 + lag_marital_status_1 +
lag_income_1 + (1 | pidp)
Data: data
REML criterion at convergence: 10242.3
Scaled residuals:
Min 1Q Median 3Q Max
-4.8643 -0.4878 0.0311 0.5199 4.7698
Random effects:
Groups Name Variance Std.Dev.
pidp (Intercept) 0.5365 0.7325
Residual 0.2623 0.5122
Number of obs: 5650, groups: pidp, 565
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.6304375 0.1501345 10.860
age 0.0156351 0.0020725 7.544
lag_retired_1yes -0.0075661 0.0633806 -0.119
lag_health_1no 0.0152044 0.0186063 0.817
lag_marital_status_1widow -0.1562918 0.0467054 -3.346
lag_marital_status_1divorce -0.1232051 0.0560884 -2.197
lag_income_1 -0.0007314 0.0031544 -0.232summ(gm2,confint = T)MODEL INFO:
Observations: 6215
Dependent Variable: PI
Type: Mixed effects linear regression
MODEL FIT:
AIC = 10260.30, BIC = 10320.05
Pseudo-R² (fixed effects) = 0.01
Pseudo-R² (total) = 0.67
FIXED EFFECTS:
------------------------------------------------------------------
Est. 2.5% 97.5% t val.
--------------------------------- ------- ------- ------- --------
(Intercept) 1.63 1.34 1.92 10.86
age 0.02 0.01 0.02 7.54
lag_retired_1yes -0.01 -0.13 0.12 -0.12
lag_health_1no 0.02 -0.02 0.05 0.82
lag_marital_status_1widow -0.16 -0.25 -0.06 -3.35
lag_marital_status_1divorce -0.12 -0.23 -0.01 -2.20
lag_income_1 -0.00 -0.01 0.01 -0.23
------------------------------------------------------------------
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
pidp (Intercept) 0.73
Residual 0.51
------------------------------------
Grouping variables:
-------------------------
Group # groups ICC
------- ---------- ------
pidp 565 0.67
-------------------------6.3 Time lag t-2
gm3 <- lmer(PI ~ age + lag_retired_2 + lag_health_2 + lag_marital_status_2 + lag_income_2 + (1 |pidp), data = data, na.action=na.exclude)
summary(gm3, correlation = FALSE)Linear mixed model fit by REML ['lmerMod']
Formula: PI ~ age + lag_retired_2 + lag_health_2 + lag_marital_status_2 +
lag_income_2 + (1 | pidp)
Data: data
REML criterion at convergence: 9231.6
Scaled residuals:
Min 1Q Median 3Q Max
-4.3741 -0.4713 0.0432 0.5158 4.6953
Random effects:
Groups Name Variance Std.Dev.
pidp (Intercept) 0.5343 0.7310
Residual 0.2564 0.5064
Number of obs: 5085, groups: pidp, 565
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.712915 0.169259 10.120
age 0.014063 0.002353 5.978
lag_retired_2yes -0.024600 0.063427 -0.388
lag_health_2no 0.027125 0.019628 1.382
lag_marital_status_2widow -0.134932 0.050824 -2.655
lag_marital_status_2divorce -0.244441 0.058991 -4.144
lag_income_2 0.005629 0.003419 1.646summ(gm3,confint = T)MODEL INFO:
Observations: 6215
Dependent Variable: PI
Type: Mixed effects linear regression
MODEL FIT:
AIC = 9249.58, BIC = 9308.38
Pseudo-R² (fixed effects) = 0.01
Pseudo-R² (total) = 0.68
FIXED EFFECTS:
------------------------------------------------------------------
Est. 2.5% 97.5% t val.
--------------------------------- ------- ------- ------- --------
(Intercept) 1.71 1.38 2.04 10.12
age 0.01 0.01 0.02 5.98
lag_retired_2yes -0.02 -0.15 0.10 -0.39
lag_health_2no 0.03 -0.01 0.07 1.38
lag_marital_status_2widow -0.13 -0.23 -0.04 -2.65
lag_marital_status_2divorce -0.24 -0.36 -0.13 -4.14
lag_income_2 0.01 -0.00 0.01 1.65
------------------------------------------------------------------
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
pidp (Intercept) 0.73
Residual 0.51
------------------------------------
Grouping variables:
-------------------------
Group # groups ICC
------- ---------- ------
pidp 565 0.68
-------------------------The value of REML is 9231.6. For the random effects, the variance of pidp intercept is 0.5343, the standard deviation is 0.7310. The variance of residual is 0.2565, the standard deviation is 0.5064. The Intra-Class Correlation ICC is 0.68. In other words, 68% of the total variability in political interest is attributed to differences between individual groups. The remaining 32% of the variability comes from differences within the same pidp group across different occasions.
The results of this thesis support the inference of hypothesis1. The estimate value of age is 0.014063, standard deviation is 0.002353. Age is statistically significant with political interest because t-value is 5.978. Holding other variable constant, an increase in age by 1 unit is associated with a change in political by 0.014063 units on average with a 95% confidence interval of (0.01, 0.02). The predicted value of political interest as figure11 showed. The relationship between age and political interest is positive. Moreover, the estimate value of intercept is 1.712915, and standard deviation is 0.169259. Intercept is statistically significant because t-value is 10.120. When all variables are 0, the average value of political interest is 1.712915 with a 95% confidence interval of (1.38, 2.04).
The estimate value of retirement is -0.024600, the standard deviation is 0.063427. Compared to those who were not retired two periods ago, those who retired two periods ago are expected to have a 0.026 units lower political interest with confidence interval (-0.15, 0.10). However, t-value is -0.388, which means retirement may not statistically significant.
The estimate value of income is 0.005629, standard deviation is 0.003419. The t-value of income is 1.646, which is statistically significant. An increase in income by 1% from two periods ago is associated with a change in political interest by 0.00005629 units on average at 90% confidence interval (-0.00, 0.01).
The estimate value of health is 0.027125, the standard deviation is 0.019628. Compared to those who were not retired two periods ago, those who retired two periods ago are expected to have a 0.027125 units lower political interest with a confidence interval (-0.01, 0.07). However, t-value is 1.382, which is not 95% statistically significant.
The estimate value of widow is -0.134932, the standard deviation is 0.050824. The t-value of widow is -2.655, which is statistically significant. Compared to other marital statuses, those who were widowed two periods ago are expected to have a 0.134932 units lower political interest, with a 95% confidence interval of (-0.23, -0.04). The estimate value of dicorce is -0.244441, the standard deviation is 0.058991, t-value is -4.141. Divorce is statistically significant. Compared to other marital statuses, those who were divorced two periods ago are expected to have a 0.244 units lower political interest, with a 95% confidence interval of (-0.36, -0.13).
stargazer(gm3, type = 'text')
=======================================================
Dependent variable:
---------------------------
PI
-------------------------------------------------------
age 0.014***
(0.002)
lag_retired_2yes -0.025
(0.063)
lag_health_2no 0.027
(0.020)
lag_marital_status_2widow -0.135***
(0.051)
lag_marital_status_2divorce -0.244***
(0.059)
lag_income_2 0.006*
(0.003)
Constant 1.713***
(0.169)
-------------------------------------------------------
Observations 6,215
Log Likelihood -4,615.788
Akaike Inf. Crit. 9,249.576
Bayesian Inf. Crit. 9,308.382
=======================================================
Note: *p<0.1; **p<0.05; ***p<0.01gm4 <- lmer(PI ~ age + lag_retired_2 + lag_health_2 + lag_marital_status_2 + log_income + (1 |pidp), data = data, na.action=na.exclude)
summary(gm4, correlation = FALSE)Linear mixed model fit by REML ['lmerMod']
Formula: PI ~ age + lag_retired_2 + lag_health_2 + lag_marital_status_2 +
log_income + (1 | pidp)
Data: data
REML criterion at convergence: 9231.8
Scaled residuals:
Min 1Q Median 3Q Max
-4.6730 -0.4711 0.0425 0.5152 4.6948
Random effects:
Groups Name Variance Std.Dev.
pidp (Intercept) 0.5343 0.7309
Residual 0.2565 0.5064
Number of obs: 5085, groups: pidp, 565
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.686239 0.171462 9.834
age 0.014449 0.002346 6.160
lag_retired_2yes -0.026014 0.063412 -0.410
lag_health_2no 0.027261 0.019631 1.389
lag_marital_status_2widow -0.135291 0.050831 -2.662
lag_marital_status_2divorce -0.244269 0.058992 -4.141
log_income 0.005621 0.003661 1.536summ(gm4,confint = T)MODEL INFO:
Observations: 6215
Dependent Variable: PI
Type: Mixed effects linear regression
MODEL FIT:
AIC = 9249.79, BIC = 9308.60
Pseudo-R² (fixed effects) = 0.01
Pseudo-R² (total) = 0.68
FIXED EFFECTS:
------------------------------------------------------------------
Est. 2.5% 97.5% t val.
--------------------------------- ------- ------- ------- --------
(Intercept) 1.69 1.35 2.02 9.83
age 0.01 0.01 0.02 6.16
lag_retired_2yes -0.03 -0.15 0.10 -0.41
lag_health_2no 0.03 -0.01 0.07 1.39
lag_marital_status_2widow -0.14 -0.23 -0.04 -2.66
lag_marital_status_2divorce -0.24 -0.36 -0.13 -4.14
log_income 0.01 -0.00 0.01 1.54
------------------------------------------------------------------
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
pidp (Intercept) 0.73
Residual 0.51
------------------------------------
Grouping variables:
-------------------------
Group # groups ICC
------- ---------- ------
pidp 565 0.68
-------------------------gm5 <- lmer(PI ~ age + lag_retired_2 + lag_health_2 + lag_marital_status_2 + lag_income_1 + (1 |pidp), data = data, na.action=na.exclude)
summary(gm5, correlation = FALSE)Linear mixed model fit by REML ['lmerMod']
Formula: PI ~ age + lag_retired_2 + lag_health_2 + lag_marital_status_2 +
lag_income_1 + (1 | pidp)
Data: data
REML criterion at convergence: 9232.3
Scaled residuals:
Min 1Q Median 3Q Max
-4.7055 -0.4681 0.0434 0.5171 4.6942
Random effects:
Groups Name Variance Std.Dev.
pidp (Intercept) 0.5346 0.7312
Residual 0.2565 0.5064
Number of obs: 5085, groups: pidp, 565
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.699077 0.170663 9.956
age 0.014314 0.002346 6.102
lag_retired_2yes -0.025706 0.063436 -0.405
lag_health_2no 0.026832 0.019628 1.367
lag_marital_status_2widow -0.135055 0.050835 -2.657
lag_marital_status_2divorce -0.243997 0.058998 -4.136
lag_income_1 0.005143 0.003859 1.333summ(gm5,confint = T)MODEL INFO:
Observations: 6215
Dependent Variable: PI
Type: Mixed effects linear regression
MODEL FIT:
AIC = 9250.27, BIC = 9309.07
Pseudo-R² (fixed effects) = 0.01
Pseudo-R² (total) = 0.68
FIXED EFFECTS:
------------------------------------------------------------------
Est. 2.5% 97.5% t val.
--------------------------------- ------- ------- ------- --------
(Intercept) 1.70 1.36 2.03 9.96
age 0.01 0.01 0.02 6.10
lag_retired_2yes -0.03 -0.15 0.10 -0.41
lag_health_2no 0.03 -0.01 0.07 1.37
lag_marital_status_2widow -0.14 -0.23 -0.04 -2.66
lag_marital_status_2divorce -0.24 -0.36 -0.13 -4.14
lag_income_1 0.01 -0.00 0.01 1.33
------------------------------------------------------------------
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
pidp (Intercept) 0.73
Residual 0.51
------------------------------------
Grouping variables:
-------------------------
Group # groups ICC
------- ---------- ------
pidp 565 0.68
-------------------------gm6 <- lmer(PI ~ age + lag_retired_1 + lag_health_2+ lag_marital_status_2 + lag_income_2 + (1 |pidp), data = data, na.action=na.exclude)
summary(gm6, correlation = FALSE)Linear mixed model fit by REML ['lmerMod']
Formula: PI ~ age + lag_retired_1 + lag_health_2 + lag_marital_status_2 +
lag_income_2 + (1 | pidp)
Data: data
REML criterion at convergence: 9231.6
Scaled residuals:
Min 1Q Median 3Q Max
-4.3729 -0.4719 0.0426 0.5161 4.6952
Random effects:
Groups Name Variance Std.Dev.
pidp (Intercept) 0.5343 0.7310
Residual 0.2564 0.5064
Number of obs: 5085, groups: pidp, 565
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.709030 0.169441 10.086
age 0.014032 0.002352 5.966
lag_retired_1yes -0.014479 0.064270 -0.225
lag_health_2no 0.027154 0.019628 1.383
lag_marital_status_2widow -0.135105 0.050826 -2.658
lag_marital_status_2divorce -0.244088 0.058993 -4.138
lag_income_2 0.005643 0.003419 1.650summ(gm6,confint = T)MODEL INFO:
Observations: 6215
Dependent Variable: PI
Type: Mixed effects linear regression
MODEL FIT:
AIC = 9249.65, BIC = 9308.46
Pseudo-R² (fixed effects) = 0.01
Pseudo-R² (total) = 0.68
FIXED EFFECTS:
------------------------------------------------------------------
Est. 2.5% 97.5% t val.
--------------------------------- ------- ------- ------- --------
(Intercept) 1.71 1.38 2.04 10.09
age 0.01 0.01 0.02 5.97
lag_retired_1yes -0.01 -0.14 0.11 -0.23
lag_health_2no 0.03 -0.01 0.07 1.38
lag_marital_status_2widow -0.14 -0.23 -0.04 -2.66
lag_marital_status_2divorce -0.24 -0.36 -0.13 -4.14
lag_income_2 0.01 -0.00 0.01 1.65
------------------------------------------------------------------
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
pidp (Intercept) 0.73
Residual 0.51
------------------------------------
Grouping variables:
-------------------------
Group # groups ICC
------- ---------- ------
pidp 565 0.68
-------------------------6.4 Compare the models
stargazer(m0, m1, m2, m3, gm1, gm2, gm3, gm4, gm5,gm6, type = 'text')
=========================================================================================================================================
Dependent variable:
-------------------------------------------------------------------------------------------------------------
PI
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
-----------------------------------------------------------------------------------------------------------------------------------------
age 0.012*** 0.012*** 0.174*** 0.013*** 0.016*** 0.014*** 0.014*** 0.014*** 0.014***
(0.002) (0.002) (0.041) (0.002) (0.002) (0.002) (0.002) (0.002) (0.002)
I(age2) -0.001***
(0.0003)
retiredyes -0.010
(0.063)
healthyes -0.004
(0.017)
marital_statuswidow -0.128***
(0.043)
marital_statusdivorce -0.050
(0.054)
log_income 0.002 0.006
(0.003) (0.004)
lag_retired_1yes -0.008 -0.014
(0.063) (0.064)
lag_health_1no 0.015
(0.019)
lag_marital_status_1widow -0.156***
(0.047)
lag_marital_status_1divorce -0.123**
(0.056)
lag_income_1 -0.001 0.005
(0.003) (0.004)
lag_retired_2yes -0.025 -0.026 -0.026
(0.063) (0.063) (0.063)
lag_health_2no 0.027 0.027 0.027 0.027
(0.020) (0.020) (0.020) (0.020)
lag_marital_status_2widow -0.135*** -0.135*** -0.135*** -0.135***
(0.051) (0.051) (0.051) (0.051)
lag_marital_status_2divorce -0.244*** -0.244*** -0.244*** -0.244***
(0.059) (0.059) (0.059) (0.059)
lag_income_2 0.006* 0.006*
(0.003) (0.003)
Constant 2.682*** 1.847*** 1.844*** -3.757*** 1.811*** 1.630*** 1.713*** 1.686*** 1.699*** 1.709***
(0.032) (0.128) (0.163) (1.426) (0.134) (0.150) (0.169) (0.171) (0.171) (0.169)
-----------------------------------------------------------------------------------------------------------------------------------------
Observations 6,215 6,215 6,215 6,215 6,215 6,215 6,215 6,215 6,215 6,215
Log Likelihood -5,550.904 -5,533.646 -5,459.713 -5,533.075 -5,543.179 -5,121.148 -4,615.788 -4,615.896 -4,616.133 -4,615.825
Akaike Inf. Crit. 11,107.810 11,075.290 10,931.430 11,076.150 11,104.360 10,260.300 9,249.576 9,249.792 9,250.267 9,249.649
Bayesian Inf. Crit. 11,128.010 11,102.230 10,971.830 11,109.820 11,164.970 10,320.050 9,308.382 9,308.598 9,309.073 9,308.456
=========================================================================================================================================
Note: *p<0.1; **p<0.05; ***p<0.01Overall, the 7th model gm3 is the best fitting model among all multilevel models.
7. Fitting, Prediction value
# Visualising the coefficients of the gm3 model
plot_model(gm3,type = 'std') +
theme_bw()
## Visualising the random effect of the gm3 model
# The figure shows the random effects of the gm3 model. Each row represents a older person
plot_model(gm3,type = 're') +
theme_bw() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank())
# Predict the value
plot_model(gm3, type = 'pred')$age
$lag_retired_2
$lag_health_2
$lag_marital_status_2
$lag_income_2